Looking into DNA breathing dynamics via quantum physics 
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We study generic aspects of bubble dynamics in DNA under time dependent perturbations, for 
example temperature change, by mapping the associated Fokker-Planck equation to a quantum 
time- dependent Schrodinger equation with imaginary time. In the static case we show that the 
eigenequation is exactly the same as that of the /3-deformed nuclear liquid drop model, without 
the issue of non-integer angular momentum. A universal breathing dynamics is demonstrated by 
using an approximate method in quantum mechanics. The calculated bubble autocorrelation func- 
tion qualitatively agrees with experimental data. Under time dependent modulations, utilizing the 
adiabatic approximation, bubble properties reveal memory effects. 
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Introduction. The stability of the double helix struc- 
ture of DNA can be attributed to the phosphodiester 
bonds in the single stranded sugar backbone and hy- 
drogen bonds between complementary base pairs of op- 
posite strands. However, the hydrogen bonds between 
parallel strands can be locally broken under physiologi- 
cal conditions preceding events such as DNA replication, 
transcription, denaturation and protein binding A 
change in environmental conditions such as pH or tem- 
perature may provide the energy required to progressively 
open the hydrogen bonds, producing domains of single- 
stranded DNA (bubbles). Eventually, e.g. upon heat- 
ing, denaturation occurs and the two strands separate 
altogether. Understanding the underlying mechanisms 
behind breathing fluctuations Q and force-assisted de- 
naturation 0], may provide further insights onto DNA 
structure and function. 

Breathing dynamics was recently detected through flu- 
orescence fluctuations in a tagged double stranded DNA 
[3| ■ Various treatments were employed for simulating this 
effect: the master-equation approach stochastic dy- 
namic simulations of the Peyrard-Bishop-Dauxois model 
@, and by adopting the Poland-Scheraga free energy 
function [7|, solving the associated Fokker-Planck equa- 
tion Q . Specifically, it has been suggested that thermally 
induced breathing processes could be mapped into the 
quantum Coulombproblem, with non-integer orbital an- 
gular momentum [9(. Here the temperature, a parameter 
in the free energy, plays a role in distinguishing repulsive 
from attractive Coulombic potentials. 

In this letter we are concerned with DNA bubble dy- 
namics when temperature, or other control parameter, 
varies in time. Adopting a generic unpairing energy func- 
tion, we study the bubble survival behavior based on 
the mapping of the Fokker-Planck equation with time- 
dependent parameters, to the quantum time- dependent 



Schrodinger equation with imaginary time. By employ- 
ing approximate quantum mechanics methods, a univer- 
sal breathing dynamics is demonstrated, insensitive to 
the details of the free energy function. Moreover, we ex- 
emplify memory effects when external parameters (e.g. 
temperature or pB) are slowly varied. 

Model. The Poland-Scheraga free energy for a single 
bubble can be written as 0, Q 



Fix) = T(x) + ck B T\n{x + 1) + 7 „, 



(1) 



with x > as the bubble size in units of base pairs. 
7o, the free energy barrier to form the initial bubble, is 
next omitted as it only introduces a constant shift in 
energy. The entropy loss associated with the formation 
of a closed polymer ring is incorporated by the factor 
ck B T\n{x + 1), whereas T(x) = 2k B T f x e(y)dy repre- 
sents the free energy for the dissociation of x base pairs 
[ToL dH , k B is the Boltzmann constant and T is the tem- 
perature. The function T(x) or s(y) may be modeled 
based on experimental data. A simple model Q assumes 
££x, where AT = T - T m with T m be- 



that r(x) = ~7l rp 

ing the melting temperature and 71 



4fc B T 7 ; T 7 = 37°C 



is the reference temperature. Since we are interested here 
in the time-evolution of the bubble distribution due to a 
change in a parameter k, e.g. temperature or pH, we 
write r = r(A/t, x) where Ak = k — k c ; k c the critical 
value of k. Note that r should be an odd function of Ak. 

At a finite temperature, the one dimensional bub- 
ble dynamics can be modeled using the overdamped 
Langevin equation with a Gaussian white noise Q 

dF 

x = -D—+ m ( V (t)n(T))=2k B TD6(t-T), (2) 

where D is a kinetic coefficient of units (k B Tx s) _1 . The 
corresponding probability density P = P(x, t) satisfies 



2 



the Fokker-Planck equation [Tj 



dt dx yj ' 2 9a; 2 ' 



(3) 



where /' = df/dx, and f(x) = F{x)/2k B T. The 

time variable was redefined 2DksTt — > t. Introducing 

-f< \~ - r(a ° 
a dressed transformation, P = e J( X >P= e 



'{x + 



l)-f*P; fi = c/2, leads to 



dP Id 2 



with a time dependent potential energy 



V(x,t) = U(x) + 



2a; 2 



a/ 
at' 



(5) 



where we assumed that the time dependent parameters 
are T and e. The potential U(x) is given by 



U(x) 



e(x) 2 jj,e{x) s'{x) 



(6) 



assuming that x ^> 1. Eq. ((4| resembles the time- 
dependent Schrodinger equation with imaginary time for 
a particle in a time-dependent potential. For the static 
case the dynamics superficially resembles the radial equa- 
tion of a particle in a central potential U (x) with cen- 
trifugal barrier + l)/2a; 2 . However, in the quantum- 
mechanical case the angular momentum fi must be an 
integer. It is thus of fundamental interest to identify a 
quantum system which permits real values for \x. 

Nuclear Liquid Drop Model. The Bohr Hamiltonian 
[l3| in the nuclear liquid drop model with a mass param- 
eter B 2 = 1 is given by (H = 1) 



i a 4 a c 5 ( 7 ,n) 
/? 4 <9/r a/3 /3 2 



+ ^09,7). (7) 



Here (3 and 7 are the parameters corresponding to 
the shape of a nucleus as an incompressible drop with 
quadrupole deformation, Q is the Euler angle onto the 
body-fixed axes, and 65(7, J7) is the Casimir operator 
of the SO(5) group [14j. For a family of potentials 
V(fi,i) = U{J3) + V(^)//3 2 [H, the /3 degree of freedom 
can be separated, 



2 a/3 2 KP) 2/3 2 
[C 5 (7,n)+27(7M7,fi) = ^ 



«(/?)= (8) 

i)-2M 7s n), 



so as the Bohr Hamiltonian eigenstates are given by 
m(/3)(^(7, fi)//3 2 . In a 7- unstable situation, ^(7) = 0, 
fi = 1,2, are integers. However, in general situations 
the effective Hamiltonian H = -§gj£r + ?7(/3) + ^^p^ 
has the exact same form as that of the breathing bubble 
with /3 replaced by x, and /x any positive number. 



This suggests that a nuclear liquid drop model, rather 
than a particle in a central potential 0] , better describes 
bubble dynamics in double-stranded polymers. 

Static limit. When all variables are time-independent 
the probability density P of can be expanded in the 
normalized eigenstates ty n solving H^ r 



P(x,t) 



(9) 



where the coefficients c n are determined by the initial 
condition and the completeness of ty n . We specify next 
the boundary conditions and distinguish between scat- 
tering potentials and binding potentials. To account 
for bubble closure absorbing boundary conditions are 
taken for vanishing bubble size, 5'„(0) = 0. Likewise, 
for considering a complete denaturation of a long strand 
with a maximum bubble size L, the absorbing condition 
^ n {L — > 00) = is implied. In order to satisfy both con- 
ditions, the family of functions e(x) should be monotonic 
for large x values so that V(x) is a binding potential. 
For instance, if e(x) is a polynomial of degree M > 0, 
the generated potential V(x) [see Eqs. ([5]) and ©] is 
always a binding potential with the asymptotic behav- 
e(x)_ x-Hx ^ ^ contrast, if M = 0, 



i°r — rrrj — — ' 00, 



00. 



2x 2 2 

e(x) is a constant corresponding to the Coulomb's poten- 
tial, and the total potential is now a scattering potential, 
allowing the function ^ n {L — > 00) to differ from zero. 

WKB Analysis. When time approaches infinity the 
transition probability ((9]) reads 



P(x 1 t)e E ^ t ^c g e-^^ g (x), 



(10) 



where ^ g (x) is the ground state of the given poten- 
tial with eigenenergy E g . In the scattering case VP g (x) 
is an oscillating function of x, while a bound ground 
state is usually nodeless and localized at a certain re- 
gion of x. What is the effect of the factor e~^ x ^ on 
the dynamics? When acting on the scattering ground 
state it affects the long time behavior of the transition 
probability leading to closure or denaturation of DNA 
bubbles 9]. On the other hand, a bound ground state 
g (x) approaches zero when x — ► 00, thus the role of 
the e~^ x ^ factor becomes influential. If the speed of 
its divergence is slower than the convergence of ^(a;), 
the bubble tends to close rather than to denaturate, and 
vice-versa. Qualitative analysis can be made in terms 
of the traditional WKB approximation The expo- 
nential factor of the ground state is given by } & g (x) oc 

e -rd Vy /2(V(y)-E g ). ^ < y( y y j n the asympt otic large 

x limit, the probability (JTTJJ) , omitting the time depen- 
dent part, reduces to P(x) oc e~ F dy [<v)+W^] 
For the Coulomb potential, e(y) = Sq is a constant, there- 
fore P oc e ~(. ea+ \ £ °\irlr[) x More generally, for bound 



potentials V(x) 



e(xY 



therefore \& g (x — * 00) 



oc 



3 



e J \ e (y}\ d y yielding the probability distribution 

P(x) ote~J* {£iv)+ ^ vWdy . (11) 

Since the integrand is non-negative, J x (s(y) + \e(y)\) dy 
either increases for s(y) > 0, leading to bubble closure, or 
does not change with x for e(y) < 0, so as the integrated 
probability linearly scales with size. The WKB analysis 
thus provides a universal long time behavior, insensitive 
to the details of the unpairing energy function. However, 
the WKB method is usually not suitable for obtaining the 
exact functional behavior, an example is provided below. 

An exactly solvable example. The transition probabil- 
ity P from an initial bubble of size xo to a bubble of final 
size x at time t is given by ^ 



P(x,x ,t) = e-^H/O 1 *) < 



-E n t 



^ n (x)y n (x ), (12) 



with the initial condition P(x, Xq, 0) = 5(x — xq). At long 
times it is approximately given by 



P(x,x Q ,t) 



12? e -/W+/(xo) e - 



E„t 



^ g {x)^ g {x Q ). (13) 



In order to simplify our analysis, we consider the fol- 
lowing expansion for the unpairing function e(x) = £ + 
2eix + 0(x 2 ). Truncating the series after the linear term 
results in T(x) — 2k B T(e x + £\X 2 ), generating the po- 



eg+2ei(2/j-l) 



see 



tential U(x) = 2efx 2 + 2e 1 e x + ^ 
Eq. ©. If £i = 0, the potential reduces to the Coulomb 
potential as in [{|. However, since the effect of E\ dom- 
inates at large distances, one should consider its contri- 
bution, for example, by using a perturbation series [l^ |. 
For simplicity we assume next that £o = 0, resulting in 
the spiked harmonic oscillator potential ([5]) 

V{x) = 2e\x 2 + £1 (2 M - 1) + ^±^, (14) 



2x 2 



with the exact ground state [l| 



* 9 (*) 



f( M +|) 



{^/2\^\xY +1 e-^ , (15) 



where E g = (3 + 2 J u)|e 1 | + £i(2/i - 1) and T(z) is the 



Gamma function. We substitute Eq. (jT5|) into (fT3|) and 
obtain 



P{x, Xq, t — > oo) 



x x- 



2(2 jgi |)m+§ (x ) 

fOi + f) 

p— (ei+|ei )ai 2 



(16) 



When £i > 0, the distribution is localized near x = 0, 
implying bubble closure. In contrast, for E\ < the 
distribution leans towards larger x values, P oc x. The 
WKB approximation (jTT|) thus produced the correct ex- 
ponential factor, but it could not provide the factor x. 



The correlation function C(t), proportional to the inte- 
grated survival probability C(t) cx f Q L P(x, Xq, t)dx\ L is 
the length of the DNA chain, can be recorded experimen- 
tally [4| . We explore next this quantity as well as the first 
passage time distribution W(t) — —dC{t)/dt. 

Results for E\ > 0. Using the sum (|12| we obtain a 
superposition of exponentially decaying functions, corre- 
sponding to various relaxation modes [4(, 



C(t) = - 



(2 £l Y+^ xl» +1 ^e +A1+ ^n + "(2ei^) 



n=0 



(17) 



Here £ = e~ 1 , and L n 2 is the associated Laguerre 
Polynomial. The first passage time distribution could 
be exactly calculated, taking the time derivative of this 
expression. At fJ, = 1/2 it has the following form 



-4eit 



exp 



1— cxp(4eii) 



[1 - exp(-4 £l £)] z 



(18) 



with VF(0) = W(oo) = 0, and a maximum in between, 
resulting in a profile similar to that obtained in [9(. The 
correlation function at = 1/2 is given by 



Cu=i(*) oc 1-exp 



,1x5 



1 — exp(4eii) 



(19) 



with the long time limit C AI= i(t) oc e 4ei *. The bubble 
lifetime is therefore given by t c — l/4ei, or r c = [4£ 1 ( / u + 
i)] _1 in general cases, see (jT7| . On the other hand, at 
short times C„_i(t) oc 1 — exp(xQ/2i). Fig. Q] presents 
the correlation function using the analytical form e(x) = 
2e\X for the unpairing free energy, and \x = 1/2, see (JJ. 
Notice that the curves at different £1, corresponding e.g. 
to different temperatures or DNA structures, follow the 
same universal temporal behavior. When presented as a 
function of a rescaled time [t — > t/ti , where C(ti) = |], 
the plots collapse into a single curve, in a good agreement 
with experiments 0] and other theoretical treatments [H, 
|2(]| . Incorporating £0 should result in a similar behavior. 

Results at ei < 0. In this case the DNA fully denatures 
at long times, and correlations diverge. At ^ = 1/2 we 
can exactly obtain the first passage time distribution 



,2 r 



8(£ia:o 
e 4|ei|t _ 1 



1 — e cx p( 4 i e i I')- 1 



(20) 



and the corresponding correlation function C(t). At long 
times both scale as L 2 . 

Time dependent effects. The adiabatic approximation 
is standardly applied to describe the dynamic of systems 
under slowly varying time dependent Hamiltonians [16 1 . 
Since the relaxation time of the bubble, order of fis [4|, 
is typically shorter than the modulation time of a pa- 
rameter k, e.g. the temperature, the quantum adiabatic 
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FIG. 1: Autocorrelation function at fi=l/2 and xq=5 [Eq. 
dSJ] for ei=0.1 (full); ei = 0.2 (dotted) and ei = 0.3 
(dashed), (inset) The curves with rescaled times t — > t/t\/2- 



approximation may be applied to describe the dynamics 
in the imaginary-time Schrodinger equation Defin- 
ing an instantaneous basis of eigenenergies H(t) \n(t)) = 

E n (t) \n(t)), we obtain (n\k) = foM*1 fc ) ; where u>kn{t) — 
Ek{t) — E n (t). In the axial representation the wave func- 
tion is written as |*(t)) = J2 n a «W e ~ & drE ^ \n[t)). 
Substituting this into the imaginary-time Schrodinger 
equation we get 



-a n {n\n) 



^ afc(t )MM e -j?- (2 i, 



Under the adiabatic approximation the coefficients a n (t) 
evolve independently from each other since couplings be- 
tween states are negligible [3] . In the present case we re- 
quire that (rt ^ fc> e -/o d ™*»M < 1. If oj kn (t) > 0, the 
exponential factor is always less than 1, while for ujk n < 
it may diverge at long times. Therefore, the applicabil- 
ity of the adiabatic approximation may be questionable 
for general instantaneous states [2l| . yet for the ground 
state it is valid as long as the standard adiabatic condi- 



tion 



(n\H\k) 



<C 1 holds. Under the adiabatic approxi- 



mation the ground state amplitude evolves according to 
hg w — a g (g\g). However, since (n\n) is zero for any one- 
dimensional real wave function, the overall function prop- 
agates as ^(x,t) ~ ^ g (x,t)e~io dTE s( T ) with V g (x,t) as 
the instantaneous solution (|15[) . Consider for example 
the potential V(x,t) = {2e\-ei)x 2 +ei(2/i-l) 



2x 2 



[see Eqs. dl])-© and (fM)) ]. which has analytical instan- 
taneous eigenstates. To simplify, we further assume that 
the system initially occupies the ground state of the po- 
tential V(x,t = 0). Under the adiabatic approximation 



P(x,t) 



(2a-2e 1 ) 



-ax 2 e -f*drE g (T) 



/2 



(22) 



with the width parameter a = E\ + {e\ — ^) 1 / 2 , and 
E g {r) = o"(r)(2// + 3) — 4ei(r). Rich information can 
be obtained due to the time-dependent evolution of s±. 



First, both the width of the distribution and the peak 
position depend on i\, the rate at which the external 
parameters (e.g. temperature) is changed. Secondly, the 
processes of increasing and decreasing the control param- 
eter may reach the same value £i, yet they may result 
in different shapes of the bubble distribution. Specifi- 

W _i_ 3 

cally, the correlation function C(t) oc (a — e\) 2 4 (1 — 
e~ CTL2 )e~ fo dTE s( T ) includes the decay factor e~ dTE si T ) 
which memorizes the different pathways that e\ (t) under- 
goes. For example, the two paths E\{i) = 1 + t/100 and 
e±(t) = 1 + i 2 /100 attain the same value at t = 1, yet 
the values of J Q drE g (r) are obviously different, yielding 
distinct characteristic decay times. We expect that this 
theoretical result could be observed experimentally. 

Summary. The dynamics of a single DNA bubble un- 
der time dependent perturbations was studied by map- 
ping the associated Fokker-Planck equation to a quan- 
tum time- dependent Schrodinger equation with imagi- 
nary time. For a generic unbinding free energy func- 
tion we analyzed bubble breathing by using the WKB 
approximation, observing a universal behavior. Specif- 
ically, a spiked harmonic oscillator potential yielded re- 
sults in qualitative agreement with experimental data. 
Under slow time dependent modulations of e.g., the tem- 
perature or pR, bubble dynamics reflects memory effects. 
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